---
output: html_document
title: "Simulation Results for Currie-MacLeod"
runhead: "Simulation Results"
author: W Bentley MacLeod
abstract: "This is an R markdown file discussion estimation and results."
date: "`r format(Sys.time(), '%B %d, %Y')`"
geometry: margin=0.5in
fontfamily: mathpazo
colorlinks: true
file: "CM-Sim-Graphs.Rmd"
---
Copyright 2019 by W. B. MacLeod, wbmacleod@gmail.com
This file is part of "CM-DoctorSim".

    "CM-DoctorSim" is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    "CM-DoctorSim" is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with "CM-DoctorSim".  If not, see <https://www.gnu.org/licenses/>.


```{r,  'Init', include = TRUE}
## Initialize R
rm(list = ls()) # clean workspace
library(rmarkdown)
library(stargazer)
library(knitr)
library(tidyverse)
library(entropy)
library(lubridate)

```
# Introduction

This document produces Tables and Figures for the analysis of the simulation data produced by
the Julia Code "RunDoctor-v2.jl".

The code reads in "CM-Sim-v2.md".

Rendering the file in Rstudio or using "render.sh" in a terminal produces a "CM-Sim-Graphics.html" with figures.
Change document type above to get pdf document.

The file "SummaryStats.md" with run parameters is also produced.
There code to time stamp this file - see below.

Individual figures are produced by un-commenting the appropriate ggplot lines below.


# Overview of the Approach

The data is generated from a Julia simulation with parameters:

## Doctor Parameters

### Doctor's Reservation Probability:

* High Reservation Probability: 0.95
* Low Reservation Probability: 0.5

### Doctor's Skill = the variance of measurement error in the doctor's assessment of the patient's condition:

* Low: 0.01
* Medium: 0.1
* High: 1

## Patient Parameters

* Degree of Complementarity in PolyPharmacy: 0.6.  If a treatment consists of two drugs, then we assume that the mean 
effect is equal to 0.6 times each single drugs effects.  Hence, we assume that a combination of two drugs is more likely
to be effective than a single drug.

* Idiosyncratic Variance in PolyPharmacy: 0.01.  If a treatment consists of two drugs, the variance of the treatment is 
assumed to be (the sum of the variances times 0.6**2) plus this term capturing idiosyncratic variance.

* Number of Patients per Doctor: 300, where all patients of each doctor are drawn from the same distribution.

## Simulation Parameters

* Months of Treatment: 36

* Number of Simulation Runs: 100

Below there are two summaries of the simulation results.  The first takes means over the all types over all treatments.  
The second shows the outcomes over time.



```{r, 'Helper Functions', echo=FALSE }
# Standard deviation function that get 5% levels

sde <- function(sdd.v) {
    n = max(length(sdd.v),2)
    qt(0.95, df = n-1) * sd(sdd.v)/sqrt(n)
}
#This function takes the decisions and computes entropy (our measure of diffusion in 
#prescribing) at a sequences of dates:
ent <- function(df.l) {
    n <- length(df.l)
    v <- df.l[[1]]
    for (i in 2:n) {
        v <- v + df.l[[i]]
    }
    m = length(v)
    return(entropy(v[1:m])/log(m))
}
```
# Overall Mean Results Over All Periods
In this section we summarize over all periods to see how physician characteristics affect mean performance.
Physicians vary in terms of skill and reservation probabilities and operate either with or without guidelines.
Here we assume that the regime without guidelines allows polypharmacy.


```{r, 'Summary Statistics', echo=FALSE }
##= Create a table of summary statistics
## Read data  
data.df <- read.csv("CM-Sim-v2.csv")
data.df$Dqual  <- factor(data.df$Dqual, levels = c('Low', 'Medium', 'High'),
                         labels = c("Low Skill", "Medium Skill", "High Skill"))
data.df$Dstyle  <- factor(data.df$Dstyle, levels = c('Low-RP', 'High-RP'), labels = c("Low Res Probs", "High Res Probs"))
data.df$Poly  <- factor(data.df$Poly, levels = c("0", "1"), labels = c("No", "Yes"))
datasumstat.df <- data.df %>% group_by(Dstyle, Dqual, Poly ) %>% summarise(
                                                                        meanE = mean(Entropy),
                                                                        meanU = mean(Utility), 
                                                                        meanO = mean(Popt),
                                                                        meanP = mean(PolyOpt),
                                                                    )

now = now()
write.csv(datasumstat.df, file = paste0("SummaryStats", ".csv")) #summary stats without time stamp.
# write.csv(datasumstat.df, file = paste0("SummaryStats", now, ".csv"))
datasumstat.df  <- datasumstat.df %>% arrange(-meanU)
datasumstat.df$meanU = datasumstat.df$meanU - min(datasumstat.df$meanU)
stargazer(as.data.frame(datasumstat.df), type = "text", summary=FALSE)

## Uncomment next line to get latex table.
## stargazer(as.data.frame(datasumstat.df), summary=FALSE, type="latex", out="CM-Summary.tex")
```
## Comments

These simulations suggest that although polypharmacy has the potential to lead to better
patient outcomes (and indeed the simulation is set up so that a combination of drugs is more 
likely to be optimal than a single drug), allowing two drug combinations does not improve 
outcomes in the shorter term.

# Performance over Time

This section explores simulations exploring how doctor skill and reservation probabilities  
affect performance over time.  Here performance is defined in terms of doctor utility, which 
is given by the number of times that the patient receives a sub-optimal treatment.
Consider first the issue of how reservation probability affects performance over time for 
physicians follow guidelines. 

Since all doctors in our model are trying to do the best they can given their skill, 
increasing entropy over time can be interpreted as a measure of improved matching between
treatments and patients.  Entropy rises over time as doctors learn about the optimal drug
choice for each patient, and patient outcomes improve.

With polypharmacy, the doctor has more possible choices, and so finding the optimum is slower.
These simulations show a case in which patient outcomes are worse in the short run in the poly-
pharmacy case (i.e. without guidelines), because of these search costs.

```{r, "Learning Over Time", echo=FALSE}
datasum.df <- data.df %>% group_by(Month, Dstyle, Dqual, Poly ) %>% summarise(
                                                                     meanE = mean(Entropy),
                                                                     meanU = mean(Utility), 
                                                                     meanO = mean(Popt),
                                                                     meanP = mean(PolyOpt),
                                                                     )

##
## Evolution over 36 months
##
sd.df  <- datasum.df %>% filter(Month < 37)

## Put performance into percentage terms - the fraction of the maximum improvement possible.

minU = min(sd.df$meanU)
sd.df <- sd.df %>% mutate(percU = 1.0 - meanU/minU)
Utime.plt <- ggplot(sd.df) +
    facet_grid(Dstyle ~ Dqual) +
    geom_smooth(aes(x=Month, y = percU, linetype = Poly), color = "black", method = "loess", formula = "y ~ x") +
    theme(legend.position = "bottom") +
    scale_y_continuous(labels = scales::percent) +
    labs(
        title = "Performance over Time",
        x = "Month",
        y = "Patient Health",
        linetype = "Does doctor violate practice guideline against poly-pharmacy? "
    )

Utime.plt
```

```{r, echo=FALSE}
## Uncomment the next line to get the pdf graph used in the paper.
## ggsave(paste0("U-56-3-years-", now, ".pdf"), plot = Utime.plt ,width = 7.5, units = "in" )

Etime.plt <- ggplot(sd.df) +
    facet_grid(Dstyle ~ Dqual) +
    geom_smooth(aes(x=Month, y = meanE, linetype = Poly), color = "black", method = "loess", formula = "y ~ x") +
    theme(legend.position = "bottom") +
    labs(
        title = "Evolution of Physician Prescription Dispersion over Time",
        x = "Month",
        y = "Physician Prescription Dispersion",
        linetype = "Does doctor violate practice guidelines against poly-pharmacy? "
    )
```

                                        # Doctor Prescription Dispersion over Time

```{r, echo=TRUE}
Etime.plt

## Uncomment the next line to get pdf of entropy graph.
## ggsave(paste0("E-56-3-years", now, ".pdf"), plot = Etime.plt ,width = 7.5, units = "in" )

```
